1. Dealing with Mortality

source("manifest.r")
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.1     ✔ readr     2.1.4
✔ forcats   1.0.0     ✔ stringr   1.5.0
✔ ggplot2   3.4.2     ✔ tibble    3.2.1
✔ lubridate 1.9.2     ✔ tidyr     1.3.0
✔ purrr     1.0.1     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

Attaching package: 'kableExtra'


The following object is masked from 'package:dplyr':

    group_rows



Attaching package: 'Matrix'


The following objects are masked from 'package:tidyr':

    expand, pack, unpack


here() starts at /Users/johngraves/Dropbox/teaching/SMDM-Europe-2023

Loading required package: forecast

Registered S3 method overwritten by 'quantmod':
  method            from
  as.zoo.data.frame zoo 

Registered S3 methods overwritten by 'demography':
  method      from 
  print.lca   e1071
  summary.lca e1071

This is demography 2.0.0.9000 



Attaching package: 'MASS'


The following object is masked from 'package:dplyr':

    select



Attaching package: 'expm'


The following object is masked from 'package:Matrix':

    expm

Learning Objectives

  • Understand options for modeling background mortality based on life tables and parametric mortality models.
  • Construct a cause-deleted life table to model cause-specific and non-cause-specific death.

Approaches to Modeling Mortality

  • Standard approach is to draw from life table data:
lt %>% 
    ungroup() %>% 
    select(-x) %>% 
    select(age,lx,dx,mx,qx) %>% 
    head(n=3) %>% 
    kable() %>% 
    kable_styling()
age lx dx mx qx
0 100000 553.22 0.00556 0.00553
1 99447 38.18 0.00038 0.00038
2 99409 23.16 0.00023 0.00023
  • lx is the total number living in a hypothetical cohort (of 100,000).
  • dx is the number of deaths in the age interval.

Approaches to Modeling Mortality

  • Standard approach is to draw from life table data:
lt %>% 
    ungroup() %>% 
    select(-x) %>% 
    select(age,lx,dx,mx,qx) %>% 
    head(n=3) %>% 
    kable() %>% 
    kable_styling()
age lx dx mx qx
0 100000 553.22 0.00556 0.00553
1 99447 38.18 0.00038 0.00038
2 99409 23.16 0.00023 0.00023
  • mx is the death rate within the interval.
  • qx is the probability of death within an interval.

Life Table Data

  • Death probabilities are useful for discrete time Markov models.
lt %>% 
    ungroup() %>% 
    select(-x) %>% 
    select(age,lx,dx,mx,qx) %>% 
    head(n=3) %>% 
    kable() %>% 
    kable_styling() %>% 
    column_spec(5, background = "#c6e6f5")
age lx dx mx qx
0 100000 553.22 0.00556 0.00553
1 99447 38.18 0.00038 0.00038
2 99409 23.16 0.00023 0.00023

Life Table Data

  • Death rates are useful for discrete event simulation (DES) models.
lt %>% 
    ungroup() %>% 
    select(-x) %>% 
    select(age,lx,dx,mx,qx) %>% 
    head(n=3) %>% 
    kable() %>% 
    kable_styling() %>% 
    column_spec(4, background = "#e8c6f5")
age lx dx mx qx
0 100000 553.22 0.00556 0.00553
1 99447 38.18 0.00038 0.00038
2 99409 23.16 0.00023 0.00023

Life Table-Based Mortality

  • A common approach is to use a life table-based lookup table in the model:
tpDn_lookup <-
    c("(34,44]" = 0.0017,  
      "(44,54]" = 0.0044,
      "(54,64]" = 0.0138,
      "(64,74]" = 0.0379,
      "(74,84]" = 0.0912,
      "(84,100]" = 0.1958)
  • Highlighted row: annual probability of death for 35-44 year old.

Life Table-Based Mortality

  • What if age bins are coarse (e.g., 0, 5, 10,… vs. 0,1,2,3,…)?
    • Assume death probability/rate is constant within age bins (and cycles)?
    • In a DES model, we need to sample a time to death.
  • Can sample from survival function in underlying life table.
  • Ages in life table often topcoded (e.g., 85, 100, etc.)

Alternative: Parametric Mortality Model

  • Fit a parametric mortality model to the life table data.
  • Reduces dimensionality to a few parameters.
  • Use these parameters to sample time of death, obtain mortality rate/probability at any arbitrary time, etc.

Obtaining Life Table Data

The R packages MortalityLaws and demography facilitate life table data aquisition and modeling.

Common sources:

Processing Life Table Data

  • Code below extracts U.S. life table data from the human mortality database (HMD).

Processing Life Table Data

  • Code below extracts U.S. life table data from the human mortality database (HMD).
hmd_usa <-
    demography::hmd.mx("USA",
                       username = "<<your user name>>",
                       password = "<<your password>>", "USA")

Processing Life Table Data

  • Code below processes HMD data and constructs a (2019) life table from it.
  • Radix (cohort size) is 100,000 individuals.

Processing Life Table Data

  • Code below processes HMD data and constructs a (2019) life table from it.
  • Radix (cohort size) is 100,000 individuals.
hmd_usa <- 
    readRDS(url("https://github.com/graveja0/SMDM-Europe-2023/raw/main/_learnr/smdm-europe-2023-cvd-model/www/usa-life-table.rds"))
mortality_year = 2019 
radix = 100000

lt = 
    hmd_usa %>% 
    demography::lifetable(.,series = "total", years = mortality_year) %>% 
    as_tibble() %>% 
    rename(age = x) %>% 
    mutate_at(vars(lx,dx), function(x) x * radix)

US Life Table (2019)

Highlighted column: Remaining life expectancy at exact age x.

lt %>% 
    select(age, lx, dx, ex, qx, mx) %>% 
    head(n=4) %>% 
    kable() %>% 
    kable_styling() %>% 
    column_spec(4, background = "#c6e6f5")
age lx dx ex qx mx
0 100000 553.221 78.963 0.00553 0.00556
1 99447 38.180 78.402 0.00038 0.00038
2 99409 23.160 77.432 0.00023 0.00023
3 99385 17.689 76.450 0.00018 0.00018

US Life Table (2019)

Highlighted column: Remaining life expectancy at exact age x.

lt %>% 
    select(age, lx, dx, ex, qx, mx) %>% 
    head(n=4) %>% 
    kable() %>% 
    kable_styling() %>% 
    column_spec(4, background = "#c6e6f5")
age lx dx ex qx mx
0 100000 553.221 78.963 0.00553 0.00556
1 99447 38.180 78.402 0.00038 0.00038
2 99409 23.160 77.432 0.00023 0.00023
3 99385 17.689 76.450 0.00018 0.00018

(We’ll try to replicate these numbers with a discrete time Markov model in a bit.)

Alive-Dead Model

Alive-Dead Model

  • As a refresher, let’s construct a very simple alive-dead model using life table death probability (qx) values.
G Alive Alive Alive->Alive 1.0-qx(t) Dead Dead Alive->Dead qx(t) Dead->Dead 1.0

Parameterize the model

params = 
  list(
    t_names = c("lifetable"),         # Strategy names. 
    n_treatments = 1,                 # Number of treatments
    s_names  = c("Alive", "Dead"),    # State names
    n_states = 2,                     # Number of states
    n_cohort = 100000,                # Cohort size
    initial_age = 0,                  # Cohort starting age    
    n_cycles = 110,                   # Number of cycles in model.  
    cycle = 1,                        # Cycle length
    life_table = lt                   # Processed HMD life table data (2019, USA)
  )

Parameterize the model

params = 
  list(
    t_names = c("lifetable"),         # Strategy names. 
    n_treatments = 1,                 # Number of treatments
    s_names  = c("Alive", "Dead"),    # State names
    n_states = 2,                     # Number of states
    n_cohort = 100000,                # Cohort size
    initial_age = 0,                  # Cohort starting age    
    n_cycles = 110,                   # Number of cycles in model.  
    cycle = 1,                        # Cycle length
    life_table = lt                   # Processed HMD life table data (2019, USA)
  )
  • No specific strategies or policies under consideration; we’re just trying to replicate mortality using life table data for a hypothetical cohort of 100,000 infants.

Parameterize the model

params = 
    list(
        t_names = c("lifetable"),         # Strategy names. 
        n_treatments = 1,                 # Number of treatments
        s_names  = c("Alive", "Dead"),    # State names
        n_states = 2,                     # Number of states
        n_cohort = 100000,                # Cohort size
        initial_age = 0,                  # Cohort starting age    
        n_cycles = 110,                   # Number of cycles in model.  
        cycle = 1,                        # Cycle length
        life_table = lt                   # Processed HMD life table data (2019, USA)
    )
  • Discrete time will be in annual cycles.

Parameterize the model

params = 
    list(
        t_names = c("lifetable"),         # Strategy names. 
        n_treatments = 1,                 # Number of treatments
        s_names  = c("Alive", "Dead"),    # State names
        n_states = 2,                     # Number of states
        n_cohort = 100000,                # Cohort size
        initial_age = 0,                  # Cohort starting age    
        n_cycles = 110,                   # Number of cycles in model.  
        cycle = 1,                        # Cycle length
        life_table = lt                   # Processed HMD life table data (2019, USA)
    )
  • The lifetable “strategy” will draw on the HMD life table data processed earlier.

Age-Dependent Transition Matrix

  • We’ll next construct a transition probability matrix using the age-specific probability of death field (qx) from the life table data.

  • To do so, we need to define a function that, for a given age and cycle length, calculates the probability of death and places this within a transition matrix.

Age-Dependent Transition Matrix

params$mP <- 
  fn_mPt(1:params$n_cycles, params)

Life Table:

lt[1,c("age","qx")]
# A tibble: 1 × 2
    age      qx
  <dbl>   <dbl>
1     0 0.00553

Transition Matrix:

Age-Dependent Transition Matrix

Life Table:

lt[1,c("age","qx")]
# A tibble: 1 × 2
    age      qx
  <dbl>   <dbl>
1     0 0.00553

Transition Matrix:

params$mP[[1]][,,"lifetable"]
       to
from    Alive   Dead     
  Alive 0.99447 0.0055322
  Dead  0       1

Age-Dependent Transition Matrix

Life Table:

lt[1,c("age","qx")]
# A tibble: 1 × 2
    age      qx
  <dbl>   <dbl>
1     0 0.00553
lt[90,c("age","qx")]
# A tibble: 1 × 2
    age    qx
  <dbl> <dbl>
1    89 0.123

Transition Matrix:

params$mP[[1]][,,"lifetable"]
       to
from    Alive   Dead     
  Alive 0.99447 0.0055322
  Dead  0       1        
params$mP[[90]][,,"lifetable"]
       to
from    Alive   Dead   
  Alive 0.87677 0.12323
  Dead  0       1

Construct the Markov Trace

sim_cohort <- function(params) {
 params$t_names %>% map(~({ 
    tr_ <- t(c("alive" = params$n_cohort, "dead" = 0))
    
    res <- do.call(rbind,lapply(params$mP, function(tp) {
        tr_ <<- tr_ %*% matrix(unlist(tp[,,.x]),nrow=params$n_state)
    }))
    res <- rbind(c(params$n_cohort,0),res) 
    dimnames(res) <- list(paste0(c(0:params$n_cycles)), params$s_names)
    res
  })) %>% 
    set_names(params$t_names)
}

Construct the Markov Trace

  • Output: named list object with full Markov trace for each state.

Construct the Markov Trace

  • Output: named list object with full Markov trace for each strategy.
trace <- 
  sim_cohort(params)
trace[["lifetable"]][1:10,]
   Alive   Dead
0 100000   0.00
1  99447 553.22
2  99409 591.40
3  99385 614.56
4  99368 632.25
5  99354 646.46
6  99340 659.67
7  99328 671.69
8  99317 682.52
9  99307 693.14

Checkpoint: Life Expectancy at Birth

Life Expectancy

  • We have a model with 110 annual cycles for a cohort of 100,000 0 year olds.
  • According to the underlying life table data, life expectancy at birth is 78.96315 years.

Life Expectancy

  • We have a model with 110 annual cycles for a cohort of 100,000 0 year olds.
  • According to the underlying life table data, life expectancy at birth is 78.96315 years.
lt %>% 
    ungroup() %>% 
    select(age,qx,ex) %>% 
    head(n=1) %>% 
    kable() %>% 
    kable_styling() %>% 
    column_spec(3, background = "#c6e6f5")
age qx ex
0 0.00553 78.963

Life Expectancy

  • Let’s verify we can replicate this with our Markov model.
lt %>% 
    ungroup() %>% 
    select(age,qx,ex) %>% 
    head(n=1) %>% 
    kable() %>% 
    kable_styling() %>% 
    column_spec(3, background = "#c6e6f5")
age qx ex
0 0.00553 78.963

Life Expectancy

  1. Define the payoff for life expectancy (1 if alive, 0 otherwise)
# Define the payoff for life expectancy
payoff_life_exp = c("Alive" = 1, "Dead" = 0)

# Calculate alive years in each cycle
life_exp_cycle = (trace[[1]] * (1 / params$n_cohort)) %*% payoff_life_exp

# Integration via alternative Simpson's method (cycle correction)
alt_simp_coef <- function(i) c(17, 59, 43, 49, rep(48, i-8),
                               49, 43, 59, 17) / 48
cycle_adj     <- function(x) sum(alt_simp_coef(length(x)) * x)

total_life_exp = cycle_adj(life_exp_cycle)

Life Expectancy

  1. Calculate life years in each cycle.1
# Define the payoff for life expectancy
payoff_life_exp = c("Alive" = 1, "Dead" = 0)

# Calculate alive years in each cycle
life_exp_cycle = (trace[[1]] * (1 / params$n_cohort)) %*% payoff_life_exp

# Integration via alternative Simpson's method (cycle correction)
alt_simp_coef <- function(i) c(17, 59, 43, 49, rep(48, i-8),
                               49, 43, 59, 17) / 48
cycle_adj     <- function(x) sum(alt_simp_coef(length(x)) * x)

total_life_exp = cycle_adj(life_exp_cycle)

Life Expectancy

  1. Define a function for cycle adjustment.1
# Define the payoff for life expectancy
payoff_life_exp = c("Alive" = 1, "Dead" = 0)

# Calculate alive years in each cycle
life_exp_cycle = (trace[[1]] * (1 / params$n_cohort)) %*% payoff_life_exp

# Integration via alternative Simpson's method (cycle correction)
alt_simp_coef <- function(i) c(17, 59, 43, 49, rep(48, i-8),
                               49, 43, 59, 17) / 48
cycle_adj     <- function(x) sum(alt_simp_coef(length(x)) * x)

total_life_exp = cycle_adj(life_exp_cycle)

Life Expectancy

  1. Sum of values (with cycle adjustment applied) is total life expectancy.
# Define the payoff for life expectancy
payoff_life_exp = c("Alive" = 1, "Dead" = 0)

# Calculate alive years in each cycle
life_exp_cycle = (trace[[1]] * (1 / params$n_cohort)) %*% payoff_life_exp

# Integration via alternative Simpson's method (cycle correction)
alt_simp_coef <- function(i) c(17, 59, 43, 49, rep(48, i-8),
                               49, 43, 59, 17) / 48
cycle_adj     <- function(x) sum(alt_simp_coef(length(x)) * x)

total_life_exp = cycle_adj(life_exp_cycle)

Life Expectancy

lt %>% 
    ungroup() %>% 
    select(age,ex) %>% 
    mutate(ex_markov = total_life_exp) %>% 
    head(n=1) %>% 
    kable(col.names = c("Age","Life Expectancy (Life Table data)",
                        "Life Expectancy (Markov-Life Table)")) %>% 
    kable_styling()
Age Life Expectancy (Life Table data) Life Expectancy (Markov-Life Table)
0 78.963 78.924
  • We can repeat this exercise for different cohort starting ages and verify things line up.

Life Expectancy

le_res <- lt %>% select(age,ex)

ages <- lt$age

le_res$ex_lt_markov = 
    ages %>% map_dbl(~({
        paramsX <- modifyList(params, list(initial_age=.x))
        paramsX$mP <- 
            fn_mPt(1:paramsX$n_cycles, paramsX)
        traceX <- 
            sim_cohort(paramsX)
        life_expX_cycle = (traceX[[1]] * (1 / paramsX$n_cohort)) %*% payoff_life_exp
        total_life_expX = cycle_adj(life_expX_cycle)
        total_life_expX
    }))

le_res %>% 
    gather(method,value,-age) %>% 
    mutate(method = factor(method, labels = c("Life Table","Markov (Life Table qx)"))) %>% 
    ggplot(aes(x = age, y = value, colour = method)) + 
    geom_point() + 
    hrbrthemes::theme_ipsum_pub(base_family = "Arial") + 
    scale_colour_aaas(name = "") + 
    theme(legend.position = "top")

Mortality Modeling

Mortality Modeling

  • Depending on context and modeling environment, you may need to model mortality parametrically.
  • The MortalityLaws package has a number of mortality models we can draw from:

Mortality Modeling

type_lut <- c(
    "1"  =  "Infant mortality" ,             
    "2"  =  "Accident hump"      ,           
    "3"  =  "Adult mortality"   ,           
    "4"  =  "Adult and/or old-age mortality",
    "5"  =  "Old-age mortality"       ,      
    "6"  =  "Full age range" 
)

laws <- 
    availableLaws()$table %>% 
    as_tibble() %>% 
    arrange(TYPE) %>% 
    mutate(TYPE = type_lut[paste0(TYPE)]) %>% 
    #filter(CODE=="HP2") %>% 
    #mutate(MODEL = gsub("\\^","||",MODEL))
    mutate(MODEL = paste0("$",MODEL,"$")) %>% 
    mutate(MODEL = ifelse(CODE == "weibull","$\\mu[x] = 1/\\sigma (x/M)^{(M/\\sigma - 1)}$",MODEL)) %>% 
    mutate(MODEL = ifelse(CODE == "invweibull","$\\mu[x] = 1/\\sigma (x/M)^{[-M/\\sigma - 1]} / [e^{(x/M)^{(-M/\\sigma)}} - 1]$",MODEL)) %>% 
    mutate(MODEL = ifelse(CODE == "perks","$\\mu[x] = [A + BC^x] / [BC^{-x} + 1 + DC^x]$",MODEL)) %>% 
    mutate(MODEL = ifelse(CODE == "wittstein","$q[x] = (1/B) A^{-[(Bx)^N]} + A^{-[(M-x)^N]}$",MODEL)) %>% 
    mutate(MODEL = ifelse(CODE == "HP","$q[x]/p[x] = A^{[(x + B)^C]} + D e^{-E \\log(x/F)^2} + G H^x$",MODEL)) %>% 
    mutate(MODEL = ifelse(CODE == "HP2","$q[x] = A^{[(x + B)^C]} + D e^{-E \\log(x/F)^2} + GH^x / [1 + GH^x]$",MODEL)) %>% 
    mutate(MODEL = ifelse(CODE == "HP3","$q[x] = A^{[(x + B)^C]} + D \\exp[-E \\log(x/F)^2] + GH^x / [1 + KGH^x]$",MODEL)) %>% 
    mutate(MODEL = ifelse(CODE == "HP3","$q[x] = A^{[(x + B)^C]} + D \\exp[-E \\log(x/F)^2] + GH^x / [1 + KGH^x]$",MODEL)) %>% 
    mutate(MODEL = ifelse(CODE == "kostaki","$q[x]/p[x] = A^{[(x+B)^C]} + D \\exp[-(E_i \\log(x/F_))^2] + G H^x$",MODEL)) %>%
    mutate(MODEL = ifelse(CODE == "opperman","$\\mu[x] = A/\\sqrt{x} - B + C\\sqrt{x}$", MODEL)) %>%
    mutate(MODEL = ifelse(CODE == "invgompertz","$\\mu[x] = [1-e^{(M-x)/\\sigma}]/[e^{(M-x)/\\sigma}-1]$",MODEL)) %>%
    mutate(MODEL = ifelse(CODE == "gompertz", "$\\mu[x] = A e^{Bx}$", MODEL)) %>%
    mutate(MODEL = ifelse(CODE == "gompertz0", "$\\mu[x] = 1/\\sigma e^{(x-M)/\\sigma}$", MODEL)) %>%
    mutate(MODEL = ifelse(CODE == "makeham", "$\\mu[x] = A e^{Bx} + C$", MODEL)) %>%
    mutate(MODEL = ifelse(CODE == "makeham0", "$\\mu[x] = 1/\\sigma e^{(x-M)/\\sigma} + C$", MODEL)) %>%
    mutate(MODEL = ifelse(CODE == "strehler_mildvan", "$\\mu[x] = K e^{-V_0 (1-Bx)/D}$", MODEL)) %>%
    mutate(MODEL = ifelse(CODE == "vandermaen", "$\\mu[x] = A+Bx+Cx^2 + I/[N-x]$", MODEL)) %>%
    mutate(MODEL = ifelse(CODE == "beard", "$\\mu[x] = e^{Bx}/[1+KA e^{Bx}]$", MODEL)) %>%
    mutate(MODEL = ifelse(CODE == "beard_makeham", "$\\mu[x] = A e^{Bx}/[1 + KA e^{Bx}]$", MODEL)) %>%
    mutate(MODEL = ifelse(CODE == "ggompertz", "$\\mu[x] = A e^{Bx}/(1+AG/B[e^{Bx}-1])$", MODEL)) %>%
    mutate(MODEL = ifelse(CODE == "vandermaen2", "$\\mu[x] = A + Bx + I/[N-x]$", MODEL)) %>%
    mutate(MODEL = ifelse(CODE == "quadratic", "$\\mu[x] = A + Bx + Cx^2$", MODEL)) %>%
    mutate(MODEL = ifelse(CODE == "kannisto", "$\\mu[x] = A e^{Bx}/[1+A e^{Bx}]$", MODEL)) %>%
    mutate(MODEL = ifelse(CODE == "kannisto_makeham", "$\\mu[x] = A e^{Bx}/[1+A e^{Bx}]+C$", MODEL)) %>%
    mutate(MODEL = ifelse(CODE == "thiele", "$\\mu[a] = A e^{-Bx} + C e^{-\\frac{1}{2} D(x-E)^2}+F e^{Gx}$", MODEL)) %>%
    mutate(MODEL = ifelse(CODE == "martinelle", "$\\mu[x] = [A e^{Bx} + C]/[1+D e^{Bx}]+K e^{Bx}$", MODEL)) %>%
    mutate(MODEL = ifelse(CODE == "rogersplanck", "$q[x] = A_0 + A_1 e^{-Ax} + A_2 e^{B(x-u)-e^{-C(x-u)}}+A_3 e^{Dx}$", MODEL))

laws1 <- laws %>% 
    filter(TYPE %in% type_lut[1:2])

laws1 %>% 
    select(-YEAR,-TYPE,-FIT,-SCALE_X) %>% 
    kable(col.names = c("Name","Model","Code")) %>% 
    pack_rows(index = table(fct_inorder(laws1$TYPE))) %>%
    cat
Name Model Code
Infant mortality
Opperman \(\mu[x] = A/\sqrt{x} - B + C\sqrt{x}\) opperman
Weibull \(\mu[x] = 1/\sigma (x/M)^{(M/\sigma - 1)}\) weibull
Accident hump
Inverse-Gompertz \(\mu[x] = [1-e^{(M-x)/\sigma}]/[e^{(M-x)/\sigma}-1]\) invgompertz
Inverse-Weibull \(\mu[x] = 1/\sigma (x/M)^{[-M/\sigma - 1]} / [e^{(x/M)^{(-M/\sigma)}} - 1]\) invweibull 

laws2 <- laws %>% 
    filter(TYPE %in% type_lut[3:4])

laws2 %>% 
    select(-YEAR,-TYPE,-FIT,-SCALE_X) %>% 
    kable(col.names = c("Name","Model","Code")) %>% 
    pack_rows(index = table(fct_inorder(laws2$TYPE))) %>%
    cat
Name Model Code
Adult mortality
Gompertz \(\mu[x] = A e^{Bx}\) gompertz
Gompertz \(\mu[x] = 1/\sigma e^{(x-M)/\sigma}\) gompertz0
Makeham \(\mu[x] = A e^{Bx} + C\) makeham
Makeham \(\mu[x] = 1/\sigma e^{(x-M)/\sigma} + C\) makeham0
Perks \(\mu[x] = [A + BC^x] / [BC^{-x} + 1 + DC^x]\) perks
Strehler-Mildvan \(\mu[x] = K e^{-V_0 (1-Bx)/D}\) strehler_mildvan
Adult and/or old-age mortality
Van der Maen \(\mu[x] = A+Bx+Cx^2 + I/[N-x]\) vandermaen
Beard \(\mu[x] = e^{Bx}/[1+KA e^{Bx}]\) beard
Beard-Makeham \(\mu[x] = A e^{Bx}/[1 + KA e^{Bx}]\) beard_makeham
Gamma-Gompertz \(\mu[x] = A e^{Bx}/(1+AG/B[e^{Bx}-1])\) ggompertz 

laws3 <- laws %>% 
    filter(TYPE %in% type_lut[5])

laws3 %>% 
    select(-YEAR,-TYPE,-FIT,-SCALE_X) %>% 
    kable(col.names = c("Name","Model","Code")) %>% 
    pack_rows(index = table(fct_inorder(laws3$TYPE))) %>%
    cat
Name Model Code
Old-age mortality
Van der Maen \(\mu[x] = A + Bx + I/[N-x]\) vandermaen2
Quadratic \(\mu[x] = A + Bx + Cx^2\) quadratic
Kannisto \(\mu[x] = A e^{Bx}/[1+A e^{Bx}]\) kannisto
Kannisto-Makeham \(\mu[x] = A e^{Bx}/[1+A e^{Bx}]+C\) kannisto_makeham 

laws4 <- laws %>% 
    filter(TYPE %in% type_lut[6]) %>% 
    filter(!(CODE %in% c("HP3","HP4","carriere2","siler")))

laws4 %>% 
    select(-YEAR,-TYPE,-FIT,-SCALE_X) %>% 
    kable(col.names = c("Name","Model","Code")) %>% 
    pack_rows(index = table(fct_inorder(laws4$TYPE))) %>%
    cat
Name Model Code
Full age range
Thiele \(\mu[a] = A e^{-Bx} + C e^{-\frac{1}{2} D(x-E)^2}+F e^{Gx}\) thiele
Wittstein \(q[x] = (1/B) A^{-[(Bx)^N]} + A^{-[(M-x)^N]}\) wittstein
Heligman-Pollard \(q[x]/p[x] = A^{[(x + B)^C]} + D e^{-E \log(x/F)^2} + G H^x\) HP
Heligman-Pollard \(q[x] = A^{[(x + B)^C]} + D e^{-E \log(x/F)^2} + GH^x / [1 + GH^x]\) HP2
Rogers-Planck \(q[x] = A_0 + A_1 e^{-Ax} + A_2 e^{B(x-u)-e^{-C(x-u)}}+A_3 e^{Dx}\) rogersplanck
Martinelle \(\mu[x] = [A e^{Bx} + C]/[1+D e^{Bx}]+K e^{Bx}\) martinelle
Carriere \(l[x] = P1 l[x](weibull) + P2 l[x](invweibull) + P3 l[x](gompertz)\) carriere1
Kostaki \(q[x]/p[x] = A^{[(x+B)^C]} + D \exp[-(E_i \log(x/F_))^2] + G H^x\) kostaki

Mortality Modeling

Key takaways:

  • MortalityLaws has a number of mortality models it can fit to life table data, depending on the context and population you are modeling.

  • You can easily experiment around with different models to see how well they fit your mortality data.

Mortality Modeling

  • Another nice feature of the package is that each mortality model type has its own defined function to calculate the mortality rate (or probability) given an age and model coefficients.

  • So rather than assume a constant mortality rate within age bins, you can model the mortality rate directly for a 40.3 year old, or a 90.0001 year old.

Mortality Modeling

Generally speaking, we need three inputs:

  • age: Ages for lifetable (e.g., 0, 1, 5, …, 100; or 0,1,2,3,…,100).
  • dx: The number of deaths between exact ages x and x+1.
  • lx: Number of survivors to exact age x.

Gompertz Model

  • Because we want to stay general (i.e., model all over the age spectrum), our first attempt will be a Gompertz model.
  • Gompertz reduces mortality to two parameters: A and B.
  • Mortalitiy rate at age \(x\):

Gompertz Model

  • Because we want to stay general (i.e., model all over the age spectrum), our first attempt will be a Gompertz model.
  • Gompertz reduces mortality to two parameters: A and B.
  • Mortalitiy rate at age \(x\):

\[ mx = A \exp(Bx) \]

Gompertz Model

ages     <- lt$age[lt$age<100]
deaths   <- lt$dx[lt$age<100]
exposure <- lt$lx[lt$age<100]

gom_fit <- MortalityLaw(
    x  = ages,      # vector with ages
    Dx  = deaths,   # vector with death counts
    Ex  = exposure, # vector containing exposures
    law = "gompertz",
    opt.method = "LF2")

Gompertz Model

# Fitted coefficients
gom_fit$coefficients
         A          B 
0.00010772 0.07534542

Gompertz Model

# Fitted coefficients
gom_fit$coefficients
         A          B 
0.00010772 0.07534542 
# Hazard rate of death for a 40 year old based on Gompertz model
# mx = A exp(Bx)
gom_fit$coefficients["A"] * exp(gom_fit$coefficients["B"] * 40 )
        A 
0.0021938

Gompertz Model

# Fitted coefficients
gom_fit$coefficients
         A          B 
0.00010772 0.07534542 
# Hazard rate of death for a 40 year old based on Gompertz model
# mx = A exp(Bx)
gom_fit$coefficients["A"] * exp(gom_fit$coefficients["B"] * 40 )
        A 
0.0021938 
# Can simply use the supplied function to calcualte 
MortalityLaws::gompertz(x = 40, par = gom_fit$coefficients)
$hx
[1] 0.0021938

$par
         A          B 
0.00010772 0.07534542 

$Sx
[1] 0.97269

Gompertz Model

# Calculate death hazard rate at age 40 from fitted Gompertz model. 
mx = gom_fit$coefficients["A"] * exp(gom_fit$coefficients["B"] * 40 )

Gompertz Model

# Calculate death hazard rate at age 40 from fitted Gompertz model. 
mx = gom_fit$coefficients["A"] * exp(gom_fit$coefficients["B"] * 40 )

# Convert mortality rate to probability
-log(1-mx)
        A 
0.0021962

Gompertz Model

# Calculate death hazard rate at age 40 from fitted Gompertz model. 
mx = gom_fit$coefficients["A"] * exp(gom_fit$coefficients["B"] * 40 )

# Convert mortality rate to probability
-log(1-mx)
        A 
0.0021962 
# Compare with life table probability (qx) at age 40:
lt[41,c("age","qx")] 
# A tibble: 1 × 2
    age      qx
  <dbl>   <dbl>
1    40 0.00208

Gompertz Model

  • Gompertz doesn’t fit all age ranges well—particularly young ages.
  • (This is a well-known fact in demography.)

Gompertz Model

plot(gom_fit)

Gompertz Model

  • If we were to focus only on adults, however, this would be a nice way to go.

Gompertz Model: 40 Year Old Cohort

ages     <- lt$age[lt$age>=40 & lt$age<100]
deaths   <- lt$dx[lt$age>=40 & lt$age<100]
exposure <- lt$lx[lt$age>=40 & lt$age<100]

gom_fit40 <- MortalityLaw(
                x  = ages,      # vector with ages
                Dx  = deaths,   # vector with death counts
                Ex  = exposure, # vector containing exposures
                law = "gompertz",
                opt.method = "LF2")

Gompertz Model: 40 Year Old Cohort

plot(gom_fit40)

Heligman-Pollard

  • While you could draw from any of the aforementioned mortality models, the Heligman-Pollard model tends to fit the entire age distribution reasonably well.

  • Easy to use: just use HP in lieu of gompertz!

Heligman-Pollard

  • While you could draw from any of the aforementioned mortality models, the Heligman-Pollard model tends to fit the entire age distribution reasonably well.

  • Easy to use: just use HP in lieu of gompertz!

ages     <- lt$age[lt$age<100]
deaths   <- lt$dx[lt$age<100]
exposure <- lt$lx[lt$age<100]

hp_fit <- MortalityLaw(
                x  = ages,      # vector with ages
                Dx  = deaths,   # vector with death counts
                Ex  = exposure, # vector containing exposures
                law = "HP",
                opt.method = "LF2")

Heligman-Pollard

plot(hp_fit)

Heligman-Pollard

Mortality rate for a 40 year old:

HP(x = 40, par = hp_fit$coefficients)
$hx
[1] 0.0020409

$par
           A            B            C            D            E           F_ 
 0.000386408  0.021271772  0.107422357  0.001025285  2.871760263 33.201213579 
           G            H 
 0.000022277  1.102611045 
x   = 40
mu1 = A^((x + B)^C) + G * H^x)
mu2 = D * exp(-E * (log(x/F_))^2))
eta = ifelse(x == 0, mu1, mu1 + mu2)
hx  = eta/(1 + eta)

Death Hazard at Age 40: Gompertz vs. HP

gompertz(x = 40, par = gom_fit$coefficients)$hx
[1] 0.0021938
HP(x = 40, par = hp_fit$coefficients)$hx
[1] 0.0020409

Alive-Dead (Modeled Mortality)

Objectives

  • Fit a parameteric mortality model to life table data.
  • Use the estimated parameters to sample death times (useful for DES modeling)
  • Use the estimated parameters as background mortality rates in discrete time Markov.

1. Sampling Death Dates for DES

  • Difficult to sample death times from life tables due to coarseness of data (e.g., five-year age bins).
  • Solution:
    1. Fit a parametric mortality model (e.g., Heligman-Pollard)
    2. Construct a survival function from modeled hazards.
    3. Sample a death time from the survival function.

1. Sampling Death Dates for DES

  • Intuition is easier if we start with a defined quantile (e.g., 25th percentile of survival)

1. Sampling Death Dates for DES

  • Intuition is easier if we start with a defined quantile (e.g., 25th percentile of survival)
u = 0.25

1. Sampling Death Dates for DES

  • We can feed the quantile (0.25) and coefficients into an inverse quantile function.
  • Straightforward to do for Gompertz model
u = 0.25
qgompertz(p = u, shape = gom_fit$coefficients["B"], rate = gom_fit$coefficients["A"])
[1] 70.467

1. Sampling Death Dates for DES

  • Let’s verify this gets us close to what we see in the life table data.
u = 0.25
qgompertz(p = u, shape = gom_fit$coefficients["B"], rate = gom_fit$coefficients["A"])
[1] 70.467
lt %>% mutate(Sx = exp(-cumsum(mx))) %>% 
  filter(age>=70 & age<=72) %>% 
  select(age,mx,Sx) %>% 
  kable() %>% 
  kable_styling() %>% 
  column_spec(3, background = "#c6e6f5")
age mx Sx
70 0.01827 0.76590
71 0.01995 0.75077
72 0.02162 0.73472

1. Sampling Death Dates for DES

  • There is no inverse quantile function for Heligman-Pollard, but we can create one!

1. Sampling Death Dates for DES

  • There is no inverse quantile function for Heligman-Pollard, but we can create one!
  • Approxfun returns a function that provides a value of y for a given value of x.
qHP <- 
  approxfun(
    y = 0:115,
    x = 1 - exp(-cumsum(HP(x = 0:115,
                           par = hp_fit$coefficients)$hx)))

1. Sampling Death Dates for DES

  • For approxfun():
    • x is quantile in CDF for survival, i.e., 1 - Survival1
    • y is the age.
  • So qHP(0.5) would return median survival age in the modeled data.

1. Sampling Death Dates for DES

qHP <- 
  approxfun(y = 0:115,
            x = 1 - exp(-cumsum(HP(x = 0:115, par = hp_fit$coefficients)$hx)))
qHP(0.25)
[1] 71.069

Survival in life table:

lt %>% mutate(Sx = exp(-cumsum(mx))) %>% 
  #filter(age>=82 & age<=83) %>% 
  filter(age>=70 & age<=72) %>% 
  select(age,mx,Sx) %>% 
  kable() %>% 
  kable_styling() %>% 
  column_spec(3, background = "#c6e6f5")
age mx Sx
70 0.01827 0.76590
71 0.01995 0.75077
72 0.02162 0.73472

1. Sampling Death Dates for DES

  • We now have the tools we need to sample death times that closely match underlying mortality data!
  • For a given simulated patient, sample a uniform (0,1) and feed this number into the inverse quantile function.

1. Sampling Death Dates for DES

  • We now have the tools we need to sample death times that closely match underlying mortality data!
  • For a given simulated patient, sample a uniform (0,1) and feed this number into the inverse quantile function.
# Sample death times for 10 individuals
u <- runif(10, min = 0, max = 1)
qHP <- 
  approxfun(y = 0:115,
  x = 1 - exp(-cumsum(HP(x = 0:115, par = hp_fit$coefficients)$hx)))
qHP(u)
 [1] 85.358 81.844 65.418 83.869 83.884 95.825 99.221 81.527 81.715 41.024

2. Modeled Mortality (Markov)

  • Let’s now incorporate modeled mortality parameters into our Alive-Dead model.
  • Rather than draw a probability of death from the life table data, we’ll calculate the probability of death at a given age using the model coefficients.

Parameterize the model

params = 
  list(
    t_names = c("lifetable",          # Strategy names
                "heligman-pollard",
                "gompertz"),          
    n_treatments = 3,                 # Number of treatments
    s_names  = c("Alive", "Dead"),    # State names
    n_states = 2,                     # Number of states
    n_cohort = 100000,                # Cohort size
    initial_age = 40,                 # Cohort starting age    
    n_cycles = 110,                   # Number of cycles in model.  
    cycle = 1,                        # Cycle length
    life_table = lt,                  # Processed HMD life table data (2019, USA)
    gom_fit = gom_fit,                # Fitted Gompertz model object
    hp_fit = hp_fit                   # Fitted Heligman-Pollard model object
  )

Age-Dependent Transition Matrices

fn_mPt <- function(t, params) {
  with(params, {
    h = cycle # Time step
    lapply(t, function(tt){
       # Life table method
        current_age_lt = pmin(initial_age  + (tt)*h - 1,max(lt$age))
        p_death <- lt[lt$age==current_age_lt,"qx"]
        # parametric mortality methods
        current_age = initial_age  + (tt)*h - 1
        r_death_gompertz <- gompertz(current_age,params$gom_fit$coefficients)$hx
        p_death_gompertz <- 1 - exp(-r_death_gompertz)
        r_death_hp <- HP(current_age, params$hp_fit$coefficients)$hx
        p_death_hp <- 1 - exp(-r_death_hp)
        
        array(data = c(1 - p_death, 0, 
                 p_death,1,
                 
                 1 - p_death_hp, 0, 
                 p_death_hp,1,
                 
                1 - p_death_gompertz, 0, 
                 p_death_gompertz,1),
              
              dim = c(n_states, n_states, n_treatments),
              dimnames = list(from = s_names,
                              to = s_names,
                              t_names)) 
    })
  })
}

Construct a Markov Trace

# Add transition matrices into params.
params$mP <- 
  fn_mPt(1:params$n_cycles, params)

# Construct the Markov trace. 
trace <- 
  sim_cohort(params)

Construct a Markov Trace

trace[["lifetable"]][1:10,]
   Alive    Dead
0 100000    0.00
1  99792  207.68
2  99584  416.43
3  99365  634.88
4  99138  861.67
5  98899 1101.19
6  98639 1361.45
7  98358 1641.68
8  98054 1945.62
9  97728 2271.99
trace[["gompertz"]][1:10,]
   Alive    Dead
0 100000    0.00
1  99781  219.14
2  99545  454.89
3  99292  708.46
4  99019  981.15
5  98726 1274.35
6  98410 1589.51
7  98072 1928.21
8  97708 2292.12
9  97317 2682.98
trace[["heligman-pollard"]][1:10,]
   Alive    Dead
0 100000    0.00
1  99796  203.88
2  99584  416.04
3  99363  637.42
4  99131  869.10
5  98888 1112.31
6  98632 1368.43
7  98361 1639.00
8  98074 1925.75
9  97769 2230.56

Calculate Life Expectancy

# Define the payoff for life expectancy
payoff_life_exp = c("Alive" = 1, "Dead" = 0)

# Calculate alive years in each cycle
life_exp_cycle = lapply(trace, function(tr) (tr*(1 / params$n_cohort)) %*%
                 payoff_life_exp)

# Integration via alternative Simpson's method (cycle correction)
alt_simp_coef <- function(i) c(17, 59, 43, 49, rep(48, i-8),
                               49, 43, 59, 17) / 48
cycle_adj     <- function(x) sum(alt_simp_coef(length(x)) * x)

total_life_exp = lapply(life_exp_cycle,cycle_adj)

Life Expectancy: Alive-Dead Model

total_life_exp
$lifetable
[1] 40.959

$`heligman-pollard`
[1] 41.008

$gompertz
[1] 41.29

Cause-Deleted Mortality

Motivation

  • Often, background mortality is drawn directly from life table data.
  • Models also frequently include cause-specific death as a separate tracked event or health state.

Motivation

  • This approach “double-counts” death because cause-specific death is also reflected in the life table estimates!
  • OK if cause-specific death is a small contributor to overall death rates.
  • But what if modeled disease is a frequent cause of death (e.g., cardiovascular disease, cancer)?

Solution

  • We can construct an alternative life table that “nets out” cause-specific deaths.
  • Not foolproof: need to make an assumption that we can parse out deaths as indepdendent contributors to overall mortality.

CVD Natural History Model

CVD Natural History Model

  • Next set of slides will walk through the process of constructing a natural history model for cardiovascular disease.
  • We’ll construct a model with non-CVD background mortality, and CVD mortality.
  • First step is constructing a cause-deleted life table.
    • Needed because at older ages, CVD is a substantial contributor (~50%) of overall mortality.

CVD-Deleted Life Table

Constructing the cause-deleted life table requires two necessary inputs:

  1. Overall mortality by age.
  2. Cause-specific mortality by age (ideally as a % of all deaths at a given age).
  • These are obtainable for many countries from the Global Burden of Disease and Human Cause of Death Database websites.

CVD-Deleted Life Table

Screenshot of Global Burden of Disease Data

We extract the percentage of overall deaths that are from CVD by age bin:

ihme_cvd <- 
  tibble::tribble(
        ~age_name,        ~val,
               1, 0.038771524,
               5, 0.038546046,
              10, 0.044403585,
              15, 0.033781126,
              20, 0.035856165,
              25, 0.053077797,
              30, 0.086001439,
              35, 0.130326551,
              40, 0.184310334,
              45,  0.21839762,
              50, 0.243705394,
              55, 0.256334637,
              60,  0.26828001,
              65, 0.272698709,
              70,  0.28529754,
              75, 0.310642009,
               0, 0.016750489,
              80, 0.353518012,
              85, 0.399856716,
              90, 0.447817792,
              95, 0.495305502
        ) %>% 
    mutate(age_ihme = cut(age_name,unique(c(0,1,seq(0,95,5),105)),right=FALSE))  %>% 
    select(age_ihme,  pct_cvd = val)

Merge these percentages into the underlying life table data and use them to calculate the total number of CVD deaths by age:

lt_ <-  # 
    lt %>% 
    mutate(age_ihme = cut(age,unique(c(0,1,seq(0,95,5),105)),right=FALSE)) %>% 
    left_join(ihme_cvd,"age_ihme") %>%
    mutate(dx_i = round(dx * pct_cvd)) %>% 
    select(age,
           age_ihme,
           pct_cvd,
           dx,  # Deaths
           dx_i, # Cause-specific deaths
           lx) %>% # Living
    mutate(a = ifelse(age==0, 0.152, 0.5)) %>% # time lived by deaths in age group
    mutate(age_interval = c(diff(age), NA)) %>% 
    select(age,age_ihme,pct_cvd, lx,dx,dx_i,a,age_interval) %>% 
    # Probability of death in interval =  deaths / total living at beginning of interval. 
    mutate(q = replace_na(1 - lead(lx) / lx, 1))  %>% 
      # Convert probability to rate. Note that we could also use q/(age_interval - q * (age_interval - a))
    mutate(m = -log(1-q)/age_interval) 
  

lt_ %>% 
  select(age,age_ihme,pct_cvd,dx,dx_i) %>% 
  filter(age %in% c(0,10,25,50,75,98))  %>% 
  kable() %>% 
  kable_styling()
age age_ihme pct_cvd dx dx_i
0 [0,1) 0.01675 553.221 9
10 [10,15) 0.04440 11.518 1
25 [25,30) 0.05308 103.706 6
50 [50,55) 0.24371 365.721 89
75 [75,80) 0.31064 1997.766 621
98 [95,105) 0.49531 1262.593 625

We next calculate the cause-deleted probability of death (qd) and death rate (md), as well as the cause-specific probability of death (qi) and death rate (mi):

Steps

  1. Cause-specific probability of death (qi): \(q \cdot dx_i / dx\)
  2. Cacluate cause-specific rates by dividing deaths of a given cause into person-years of exposure.
  • This is equivalent to multiplying the overall rate by the ratio of deaths of a given cause to the total.

Steps

  1. Cause-specific probability of death (qi): \(q \cdot dx_i / dx\)
  2. Cacluate cause-specific rates by dividing deaths of a given cause into person-years of exposure.
    • This is equivalent to multiplying the overall rate by the ratio of deaths of a given cause to the total.
lt_ <- 
  lt_ %>% 
    # Cause specific probability of death. 
    mutate(qi = q * dx_i / dx) %>% 
    mutate(Rd = (dx - dx_i) / dx) %>% 
    mutate(md = m * Rd) %>% 
    mutate(mi = m - md) %>% 
    mutate(dxd = dx - dx_i)

lt_ %>% 
  select(age, dx, dxd, dx_i, q, qi, md, mi) %>% 
  head() %>% 
  kable() %>% 
  kable_styling()
age dx dxd dx_i q qi md mi
0 553.221 544.221 9 0.00553 0.00009 0.00546 0.00009
1 38.180 37.180 1 0.00038 0.00001 0.00037 0.00001
2 23.160 22.160 1 0.00023 0.00001 0.00022 0.00001
3 17.689 16.689 1 0.00018 0.00001 0.00017 0.00001
4 14.209 13.209 1 0.00014 0.00001 0.00013 0.00001
5 13.213 12.213 1 0.00013 0.00001 0.00012 0.00001

Modeling CVD-Deleted Mortality

  • We can also apply the methods covered earlier to the cause-deleted life table data!
ages_     <- lt_$age[lt_$age<100 & lt_$age>=25]
deaths_   <- lt_$dxd[lt_$age<100 & lt_$age>=25] 
exposure_  <- lt_$lx[lt_$age<100 & lt_$age>=25]

mort_fit_CVDdeleted <- MortalityLaw(
    x  = ages_,
    Dx  = deaths_,   # vector with death counts
    Ex  = exposure_, # vector containing exposures
    law = "HP2",
    opt.method = "LF2")

params = 
  list(
    t_names = c("natural_history"),       # Strategy names. 
    n_treatments = 1,                     # Number of treatments
    s_names  = c("Healthy", "CVD","CVDDeath","Dead"),
    n_states = 4,                         # Number of states
    n_cohort = 100000,                    # Cohort size
    initial_age = 25,                      # Cohort starting age    
    n_cycles = 100,                       # Number of cycles in model.  
    cycle = 1,                            # Cycle length
    r_H_CVD = 0.1,                        # Healthy -> CVD transition rate
    u_H = 1,                              # Utility state payoffs
    u_CVD = 1,
    u_D = 0,
    u_DCVD = 0,
    c_H = 0,                              # Cost state payoffs
    c_CVD = 0,
    c_DCVD = 0,
    c_D = 0, 
    cause_specific_mortality =            # Cause specific morality rate
      approxfun(lt_$age, lt_$mi,rule = 2),
    mort_fit = mort_fit_CVDdeleted        # Cause-Deleted Mortality Fit 
  )

Basic CVD Model

G Healthy Healthy Healthy->Healthy CVD CVD Healthy->CVD r_HS Dead Dead Healthy->Dead r_HD(t) CVD->CVD Dead(CVD) Dead (CVD) CVD->Dead(CVD)    r_DCVD(t)         CVD->Dead  r_HD(t) Dead(CVD)->Dead(CVD) Dead->Dead

Basic CVD Model

  • Healthy Death based on modeled (Heligman-Pollard) cause-deleted mortality.
G Healthy Healthy Healthy->Healthy CVD CVD Healthy->CVD r_HS Dead Dead Healthy->Dead r_HD(t) CVD->CVD Dead(CVD) Dead (CVD) CVD->Dead(CVD)    r_DCVD(t)         CVD->Dead  r_HD(t) Dead(CVD)->Dead(CVD) Dead->Dead

Basic CVD Model

  • CVD Death based on age- and cause-specific death rates from cause-deleted life table.
G Healthy Healthy Healthy->Healthy CVD CVD Healthy->CVD r_HS Dead Dead Healthy->Dead r_HD(t) CVD->CVD Dead(CVD) Dead (CVD) CVD->Dead(CVD)    r_DCVD(t)         CVD->Dead  r_HD(t) Dead(CVD)->Dead(CVD) Dead->Dead

Basic CVD Model

  • We also allow for a CVD Death transition based on non-CVD causes.
G Healthy Healthy Healthy->Healthy CVD CVD Healthy->CVD r_HS Dead Dead Healthy->Dead r_HD(t) CVD->CVD Dead(CVD) Dead (CVD) CVD->Dead(CVD)    r_DCVD(t)         CVD->Dead  r_HD(t) Dead(CVD)->Dead(CVD) Dead->Dead

Age-Dependent Transition Matrices

fn_mPt_CVD <- function(t, params) {
  with(params, {
    h = cycle # Time step
    lapply(t, function(tt){
       # Life table method
        current_age = initial_age  + (tt)*h - 1
        
        r_H_D <- HP(current_age, params$mort_fit$coefficients)$hx
        r_CVD_D = params$cause_specific_mortality(current_age) 
        
        R <- array(data = c(
                      -(r_H_D+params$r_H_CVD),0,0,0, 
                      params$r_H_CVD,-(r_CVD_D+r_H_D),0,0, 
                      0,r_CVD_D,0,0,
                      r_H_D,r_H_D,0,0),
                      
              
              dim = c(n_states, n_states, n_treatments),
              dimnames = list(from = s_names,
                              to = s_names,
                              t_names)) 
        P = apply(R,3,expm, simplify=FALSE)
        P
    })
  })
}
params$mP <- fn_mPt_CVD(1:params$n_cycles,params)

Construct a Markov Trace

sim_cohort <- function(params) {
  
 params$t_names %>% map(~({ 
    tr_ <- t(c(params$n_cohort,rep(0,params$n_states - 1)))
    
    res <- do.call(rbind,lapply(params$mP, function(tp) {
        tr_ <<- tr_ %*% matrix(unlist(tp[[.x]]),nrow=params$n_state)
    }))
    res <- rbind(c(params$n_cohort,rep(0,params$n_states - 1)),res) 
    dimnames(res) <- list(paste0(c(0:params$n_cycles)), params$s_names)
    res
  })) %>% 
    set_names(params$t_names)
}
trace <- 
  sim_cohort(params)
trace[["natural_history"]][1:10,]
  Healthy     CVD CVDDeath    Dead
0  100000     0.0   0.0000    0.00
1   90391  9506.2   0.2945  103.00
2   81702 18088.0   1.1404  208.58
3   73847 25833.6   2.4853  316.88
4   66745 32822.6   4.2817  428.11
5   60324 39126.6   6.8544  542.45
6   54519 44809.4  11.5751  660.14
7   49270 49930.7  17.3944  781.42
8   44526 54544.1  23.8195  906.57
9   40236 58697.0  31.3741 1035.87

Good News

  • Modeled CVD and non-CVD mortality closely matches the cause-deleted life table data!
lt_compare <- 
    lt_ %>% 
    mutate(dx = dx - dx_i) %>% 
    select(age = age, dead_lt = dx, cvddeath_lt = dx_i) %>% 
    mutate(age = age + 1) %>% 
    filter(age>=params$initial_age) %>% 
    mutate(dead_lt = ifelse(age==25,0,dead_lt),
           cvddeath_lt = ifelse(age==25,0,cvddeath_lt)) %>% 
    mutate_at(vars(dead_lt,cvddeath_lt), cumsum)


trace[["natural_history"]] %>% 
    data.frame() %>% 
    tibble() %>% 
    set_names(c("healthy","cvd","cvddeath","dead")) %>% 
    mutate(age = params$initial_age  + (row_number())*params$cycle - 1 ) %>% 
    select(age,cvddeath,dead) %>% 
    inner_join(lt_compare, "age") %>% 
    gather(state,value,-age) %>% 
    mutate(state = factor(state,levels = c("cvddeath","cvddeath_lt","dead","dead_lt"), 
                          labels = c("CVD Deaths",
                                     "CVD Deaths (Life Table)",
                                     "Non-CVD Mortality",
                                     "Non-CVD Mortality (Life Table)"))) %>% 
    mutate(state1 = ifelse(grepl("CVD Deaths",state),"CVD Deaths","Non-CVD Mortality")) %>% 
    mutate(state2 = ifelse(grepl("Life Table",state),"Life Table","Markov")) %>% 
    ggplot(aes(x = age, y = value, colour = state2)) + 
    geom_line(aes(lty=state1))  +
    hrbrthemes::theme_ipsum(base_family = "Arial") + 
    scale_colour_aaas(name="") + 
    scale_linetype(guide="none") +
    geom_dl(method = list("last.bumpup",hjust=1),aes(label=state1),colour = "black") +
    theme(legend.position="top") + 
    labs(x = "Age", y = "Deaths") + 
    scale_y_continuous(limits = c(0,radix))

Your Turn

Case Study Overview

TK

Learning Objectives

TK

Back to Website